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Abstract. 

Enzymatic molecules that actively support many cellular processes, including 
transport, cell division and cell motility, are known as motor proteins or molecular 
motors. Experimental studies indicate that they interact with each other and they 
frequently work together in large groups. To understand the mechanisms of collective 
behavior of motor proteins we study the effect of interactions in the transport 
of molecular motors along linear filaments. It is done by analyzing a recently 
introduced class of totally asymmetric exclusion processes that takes into account the 
intermolecular interactions via thermodynamically consistent approach. We develop a 
new theoretical method that allows us to compute analytically all dynamic properties 
of the system. Our analysis shows that correlations play important role in dynamics of 
interacting molecular motors. Surprisingly, we find that the correlations for repulsive 
interactions are weaker and more short-range than the correlations for the attractive 
interactions. In addition, it is shown that symmetry of interactions affect dynamic 
properties of molecular motors. The implications of these findings for motor proteins 
transport are discussed. Our theoretical predictions are tested by extensive Monte 
Carlo computer simulations. 
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1. Introduction 


Motor proteins or molecular motors are enzymatic molecules that actively participate 
in all major biological processes such as cellular transport, cell division, transfer of 
genetic information, synthesis of proteins, cell motility and signaling [l]-j6j. They 
transform chemical energy from specific reactions that they catalyze (usually, hydrolysis 
or biopolymerization) into mechanical work to support their functions. For example, the 
directed motion along linear cytoskeleton filaments by kinesin, myosin and dynein motor 
proteins is fueled by the hydrolysis of adenosine triphosphate (ATP) |6|. Biological 
molecular motors have been intensively studied in recent years, and currently the single- 
molecule dynamics of motor proteins is well described (4|[5||7|. Although the properties 
of individual molecules are very useful, in biological systems motor proteins typically 
function in large teams. This underlines the importance of understanding the collective 
behavior of molecular motors j5j;8j[9|. 

Experimental studies of kinesin motor proteins moving along microtubules indicate 
that these molecular motors interact with each 10-12 . It was argued that these 


interactions most probably are short-range and relatively weak attractive (1.6±0.5fc^T) 
10 . It is reasonable to assume that many other motor proteins have similar 
properties. At microscopic level, molecular motors are involved in a variety of 
chemical transitions such as binding to the filament, chemical transformations during 
the hydrolysis, dissociation from the track |6|. Intermolecular interactions influence all 
these processes, suggesting an important role for interactions in the collective behavior of 
molecular motors. However, the underlying mechanisms are still not well clarified [5jj8]. 
Existing theoretical studies of cooperative dynamics in interacting molecular motors are 
mostly phenomenological without quantitative description of relevant chemical processes 


13-16 


Recently, we introduced a new theoretical approach for analyzing collective 
properties of interacting molecular motors jlT]. It is based on using a class of non¬ 
equilibrium models called totally asymmetric simple exclusion processes (TASEP), 
which are very powerful for studying multi-particle dynamic phenomena 18 201. There 
are many processes in Chemistry, Physics and Biology that have been successfully 
analyzed by utilizing the asymmetric exclusion processes 19 -25 . TASEPs were also 


employed before for investigating dynamic properties of motor proteins |9,14,16,191,22, 
26 , including interacting molecular motors [I4| , [l6 ,27,28|. But the main advantage of 
our method is a procedure that describes all chemical transitions at the single-molecule 
level using fundamental thermodynamic concepts 17 . This allows us to properly couple 


microscopic properties of interacting molecular motors with their collective dynamic 
behavior. 

Analyzing this theoretical approach, it was found that there is an optimal 
interaction strength (weakly repulsive) that can maximize the current through the 
system |17|. In addition, the calculations suggested that correlations play important role 
in dynamics of interacting molecular motors. However, the progress in understanding the 
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cooperativity in motor proteins in this model was limited by the following issues. Two 
mean-field analytical treatments were proposed. But the first one, a simple mean-field 
approach (SMF), failed completely, as compared with extensive Monte Carlo computer 
simulations, producing unphysical trends in dynamic properties for strong interactions. 
The second approach, a cluster mean-field (CMF), worked better, but it involved very 
heavy numerical calculations. At the end, CMF was able to reproduce only qualitatively 
dynamics of interacting molecular motors and not even for all ranges of parameters. 
Furthermore, only symmetric splitting of interactions on hopping rates was considered. 
In addition, it was not possible to extend CMF to take into account more realistic 
features of motor protein’s transport such as backward steps, bindings and unbindings 
from the filament, and more general symmetries of the interactions. To understand 
the mechanisms of cooperativity, it is important to have an analytical method that can 
successfully capture main features of interacting molecular motors, and which can be 
also extended to more complex situations. 

In this paper, a new theoretical framework for analyzing complex dynamics of 
interacting molecular motors via TASEP is presented. We develop a modified cluster 
mean-field (MCMF) approach that accounts for some correlations in the system. This 
provides a direct way of analytically calculating all dynamic properties in the system, 
and the results agree quite well with computer simulations. The method allows us to 
explicitly analyze the role of interactions in dynamics of interacting molecular motors. 
More specifically, it is found that correlations are weaker and more short-range for 
repulsive interactions while for attractions they are stronger and more long-range. 
We also investigate the role of the symmetry of interactions and show that it might 
dramatically modify the dynamic behavior. But most importantly, the developed 
framework allows us to understand the microscopic origin of dynamic phenomena in 
motor proteins and it can be easily generalized to account for more complex processes 
associated with molecular motors. 

2. Theoretical Description 

2.1. Model 

In our model, the transport of molecular motors along linear filaments is viewed as a 
motion of multiple particles on a lattice with L (L > 1) sites, as shown in Fig. 1. The 
state of occupancy for each lattice site i (1 < i < L) is characterized by an occupation 
number t*. If the site i is occupied then t* = 1, if it is empty then r* = 0. Each lattice 
site can accommodate only one particle. 

In addition to exclusions, molecular motors can interact with each other via a short- 
range potential. Here we assume that any two neighboring particles interact with each 
other with an energy E. The case of positive E defines attractive interactions, while 
E < 0 corresponds to repulsions. In other words, any bond connecting two neighboring 
particles on the lattice is associated with the energy E. In our system transition rates 
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depend if these bonds are broken or created. Any forward motion of the individual 
molecular motor that does not change the number of bonds is taking place with the 
rate 1. It can be done by a single molecule that do not have any neighbors (see Fig. 1), 
or it can involve breaking one bond and creating another one (Fig. 1). In both cases, 
there is no energy change in the system. However, the forward transition associated 
with creating a new bond has a rate q 1. In this case, the molecular motor joins an 
existing cluster of particles: see Fig. 1. Similarly, the transition that is coupled with 
breaking the bond has a rate r 1. Here the particle dissociates from the cluster but 
simultaneously it does not bind to another cluster ( Fig. 1). The transition rates q and 
r are associated with changes in energy. 

It has been argued that creating and breaking such bonds (or pairs of particles) can 
be viewed as opposing chemical transitions, which justifies the application of the detailed 
balance arguments ED This leads to the following relation between the transitions 
rates, 


Q 

- = exp 
r 



( 1 ) 


To evaluate dynamic properties of molecular motors we need to know the explicit values 
for rates q and r. The interaction energy can be split in the following way, 


q = exp 



r = exp 


((6-l)E\ 
V k B T )’ 


( 2 ) 


where a dimensionless parameter 6 (0 < 9 < 1) specifies how the energy affects these 
transition rates. Previously, only a symmetric splitting of interactions (6 = 1/2) has 


been considered 17 


It is easy to understand the physical meaning of Eq. (|2]) [l7|. When the interactions 
between molecular motors are attractive (E > 0), the rate of creating the bond is larger 
(q > 1), while the rate of breaking the bond is smaller (r < 1). For repulsive interactions 
(E < 0) the trend is opposite — it is faster to break the particle cluster (r > 1) than to 
increase the cluster size (q < 1). In the case of no interactions (E = 0) these transitions 
rates are the same (q — r — 1) and the model reduces to standard TASEP with only 
hard-core exclusions. 


P 



binding unbinding 


Figure 1 . Schematic view of the TASEP model for interacting molecular motors. 
Binding corresponds to particle joining the cluster, while unbinding describes the 
breaking from the cluster. 
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In our model, particles enter the system from the left side of the lattice and they 
leave the system from the right side of the lattice: see Fig. 1. Interactions are also 
modifying the entrance and exit rates in comparison with the original TASEP model. 

When entering the system does not lead to creating a pair of particles the rate for this 
process is a (Fig. 1). However, entering with creating the bond has a rate qa. Similarly, 
the exit rate for the case when no bond breaking is involved is equal to /3, while the 
dissociating from the cluster is taking place with the rate r/3 (Fig. 1). 


2.2. Modified Cluster Mean-Field Theory 


Previous theoretical studies indicated that correlations are important for the system 
with multiple particles 17 . Neglecting correlations leads to unphysical behavior for 
strong interactions between molecular motors 17 . This indicates that any successful 


theoretical treatment must take correlations into account. This is the main idea of our 
approach that we call a modified cluster mean-field. 

To account for correlations we analyze bulk clusters of two neighboring sites on 
the lattice. Each cluster can be found in one of four possible states. We label a 
configuration with 2 empty sites as (0,0), with two occupied sites as (1,1), and two 
half-occupied clusters are labeled as (1,0) and (0,1). Next, we introduce functions Pu, 
Pio, Poi and Poo as probabilities of finding the configurations (1,1), (1,0), (0,1) and 
(0,0), respectively. The conservation of probability for these functions requires that 

P\\ + Pro + Poi + Poo — 1- (3) 

In addition, we have 


Pio + Pu — Pi Poi + Pn — Pi (4) 

where p is the bulk density (or the probability to find the particle at given site). Here 
it was assumed also that the bulk density is uniform and independent of the position 
on the lattice if the two-site cluster is far away from the boundaries. Combining Eqs. 
@ and @ leads to Pio = Poi and 

Pio + Poo + P — 1- (5) 


Let us consider a particle flux in the bulk of the system at large times when 
stationary conditions are achieved. We can concentrate on the segment of 4 consecutive 
sites as shown in Fig. 2. To measure the current, only transitions between the second and 
the third sites of the segment are counted. Then there are four possible configurations 
that support the particle current, as illustrated in Fig. 2. Correspondingly, the total 
flux have 4 contributions from each configurations, fl u ik = J\ + J 2 + J 3 + -h- The first 
contribution from configuration 1 (Fig. 2) can be written as 


■h = 7P10 


Pqo 

p +Poo 



( 6 ) 


where 7 = - \ >. . This expression is an approximation and it can be understood 

1+ex P{ktfr) 

in the following way. The second factor (Pio) gives the probability that the cluster 
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1 



Figure 2. Four-sites bulk lattice segments that are utilized for calculating the particle 
currents in the system. 


consisting of the second and third sites (see Fig. 2) is in the state (1, 0). The first factor 
(7) gives a probability that the first site is empty, i.e., it is just a Boltzmann’s factor. 
The third factor ( 77 ^ 7 ) is a probability to have the last site empty. If we have the 
configuration ( 1 , 0 ) in the middle cluster then the last site can be found in one of two 
states: it can be occupied with the probability p or it can be empty with the probability 
Poo- This is because in this case the cluster consisting of the sites 3 and 4 can only 
be found in configurations (0,0) or (0,1). Then the particle current from the second 
configuration (Fig. 2) is equal to 

-Poo 


J 2 = (1 - i flP 


10 


(7) 


Here (1 — 7) = 


ex P l fcfr 


P + -Poo 

is a probability to have the first site occupied, and r 


i+exp \jffff 

is the transition rate for this configuration. Similar arguments can be presented for 
contributions from configurations 3 and 4, yielding 


•h = 7<?Pio 


P 


P + Pc 


00 


and 


Ji = { 1 - 7 )^ 


10 


P 


( 8 ) 


(9) 


JJ + Poo 

Combining together Eqs. ®. 0 .® and ([9]) we obtain the expression for the total 
bulk current, 

J ( -P10-P00 \ n X f -P10-P00 

Jbulk = 7 I , n ) + (1 - l) r 


P + Poo / 
P-P 10 


P + P( 


00 


P + Poo 
P-P 10 


P + Pc 


00 


+ <?7 


+ (1 - 7 ) 


( 10 ) 
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This equation can be also written in the following form, 

P10P00 \ , r, f pp 10 


Ibulk A 


1 -P 


10 


+ B 


1-P, 


10 


where auxiliary functions A and B are defined as 

q + exp 


A = 


1 + rex P ( kCf 


1 + exp 


\k B T ) 


B = 


E 

k B T 


1 + exp 


(-E-') 

\k B T ) 


( 11 ) 


( 12 ) 


To calculate explicitly dynamic properties in the system we have to express 
everything in terms of the bulk density p and the interaction energy E. Eq. (JH]) gives 
the connection between P w , P 00 and p, and one more additional relation is needed in 
order to have all equations only in terms of p and E. We can approximate the function 
Pro as 

P(1 ~P) 


Pi o - 


p + pexp(^) 


(13) 


The physical meaning of this approximation can be explained if we note that Pio is the 
probability to have the two-site cluster in the configuration (1,0). This probability is 
equal to the product of two terms: one is the probability to have the first site occupied 
(p) and the second term is the probability that the second site is empty given that the 
first one is not ((1 — p)/[l — p + pexp (j]). It can be argued that the situation when 
two sites in the cluster are occupied is affected by the interactions between them, and we 
approximate it via the usual Boltzmann’s factor. One can see that this equation leads to 
a very reasonable behavior at the limiting cases. When there is no interactions (E = 0) 
it predicts that P w = p( 1 — p), as expected. For very strong repulsions (E —> —oo) it 
gives Pio = p, which is a correct result since in this limit our problem is identical to 
the motion of non-interacting dimers on the lattice IT, 30 . For very large attractions 


(E —> oo) the prediction is that Pio —>■ 0. This is again seems to be a reasonable result 
because in this case one is expecting to have the whole system fully occupied without 
any vacancies. 

Finally, taking into account all approximations we obtain the general expression for 
the bulk current only in terms of the particle density p and the interaction E, 


Jbulk r 


Ap{ 1 - p) 2 

1 —2P + Pexp(j f ) 


(1 -P) 2 + pexp ( 

k BT )_ 

1 - p + p exp 

?r| 

i_i 


+ 


Pp 2 (l-p) 


TT ( 14 ) 

(1 -p) 2 + pexp (^J 

For the case of zero interactions this equation suggests that fluik = p(l ~ p), the known 
result form the original TASEP model 
predict that 

P( 1 - 2p) 


18,19 . For strong repulsions {E —>■ — oo) we 


Jbulk 


P 


(15) 
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This is identical to the expression that was derived earlier for TASEP of dimers 30 


For large attractions (E —> oo) it predicts that the bulk current vanishes, .fluik = 0, 
and this is expected because particles will not be able to move since they will be stuck 
together in one large cluster. 

At boundaries the dynamics in the system is governed by exit and entrance rates. 
Using the same approximations as explained above for the bulk fluxes it can be shown 
that the expressions for entrance current is given by 


'Ien.tr 


a(l~p) 

1-2P + Pexp(^) 

+ aqp{ 1 - p) 

l-P + Pex P (jr 



(16) 


For E — 0 this equation reduces to J entr = a(l — p), which is expected for this situation. 
For E —»■ — ooit predicts that J entr = a(l — 2 p), in agreement with known results 
on TASEP of extended objects 30 . For strong attractions (E —y oo) the current 


disappears, J entr —> 0 . Similarly, for the exit current we obtain 

Pp 


•Iexit. 


1 -p + rpexp ( j-f 


1 - p + pexp 


\k B T ) 


(17) 


Again, for E = 0 and for E —> —oo it produces the expected result, J exit = j3p , while 
for strong attractions it leads to J ex it —> 0 . 


2.3. Phase Diagrams 

Similarly to the original TASEP, it can be argued that in the system of interacting 
molecular motors there are three dynamic phases at stationary conditions. When the 
rate limiting step is the entrance into the system we have a low-density (LD) phase. 
For the case of exiting being small a high-density (HD) phase will be realized. Finally, 
when bulk processes are the most important, the system is in a maximal-current (MC) 
phase. 

Our analytical theory can calculate explicitly the phase boundaries. The MC phase 
is characterized by a condition that 9J g" lk = 0 , which leads to the following expression: 

(2p 2 — 4p+ l) (p- l) 4 -p 4 exp ^(0 + 3 )~j^j -P 2 (2p- l)(p- l) 2 exp (^0 + 2)^^J 
+ (P - l) 6 exp - (P - 2)p(p - l) 4 exp ((0 + 1)^ + (p 4 - 2p 5 ) exp (4 

~P (V - 16p 2 + 15p - 4) (p - l) 2 exp - p 3 (p 3 - 10p 2 + 13p - 4) exp 

+p 2 (3p 4 - 22p 3 + 39p 2 - 26p + 6 ) exp = 0. (18) 

For E = 0 this complicated expression reduces to the following formula, 

- 4p 5 + 10p 4 - 16p 3 + 14p 2 - 8p + 2 = 0, 


( 19 ) 
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which has only one real root, p — For very large repulsions (E — > —oo) one can 
obtain from Eq.([l8|), 


(2p 2 — 4p + l) (p — l ) 4 = 0. (20) 

This equation has three roots but only one of them is physically reasonable, p = 1—1/ V2, 
which leads to a nonzero flux in the system. (The root p = 1 does not support the 
nonzero flux through the system and it can be neglected.) Substituting this density 
into Eq. (15) leads to a prediction that the particle flux in this case is equal to 
J = 3 — 2\/2 ~ 0.17. For large attractions (E — > oo) Eq.(18) predicts that p = 1/2, but 
the current here is approaching zero, as was discussed above. For general conditions, 
this equation can be always solved numerically and after choosing the physically relevant 
root for the density in the MC phase, Pmc, the particle fluxes can be calculated using 
Eq. (14). The phase diagrams calculated via this method are presented in Fig. 3 for 
various sets of parameters. 





Figure 3. Stationary phase diagrams of TASEP with interacting particles for various 
interaction strengths and interaction splittings: (a) E = — 1.7k B T, 9 = 0.25; (b) 
E = -1.2k B T, 9 = 0.5; (c) E = -0.7k B T, 6 = 0.75; (d) E = 3k B T, 9 = 0.5. Lines are 
theoretical predictions, symbols are from Monte Carlo computer simulations. 


To determine the density of molecular motors in the LD phase and the boundary 
lines separating the low-density and the maximal-current phases, we use the continuity 
of the stationary currents at the transition line, Jbuik = Jentr■ Combining Eqs. (14) and 
(16) yields the following expression for the entrance rate a, 


Ap( 1 - p) 

1 — 2 p + pex p 

Cht) 

+ Bp 2 

1 - p + p exp ^ 

k B T) 



1 — 2p + pexp(^ 

+ qp 

_(1-P) 2 + Pexp(jrfy) 




a = 


( 21 ) 
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Solving this equation for p provides an estimate for the particle density in the bulk of 
the system in the LD phase. Increasing the entrance rate a leads to large bulk densities, 
and the phase boundary between LD and MC phase is achieved when p = pmc- For 
example, for zero interactions, E = 0, from Eq. (21) it follows that pld = ot, and the 
phase boundary corresponds to a = 0.5. These estimates fully agree with results from 
the original TASEP with non-interacting particles 18,19 . For the case of very strong 
repulsions, E —» — oo, we derive from Eq. (21) that Pld = «/(1 + a), and the phase 
boundary between LD and MC phase corresponds to a = — 1 « 0.41. This again 

agrees with known results on TASEP of extended objects |30|. In the opposite limit 


of very large attractions, E —> oo, the low-density phase cannot be realized for any 
nonzero values of a. 

Similar calculations can be performed for obtaining properties of the HD phase and 
the boundaries between high-density and maximal-current phase. The exit rate f3 is 
coupled with the density p via 


p = 


A{l-pf 

i—2 P + p exp 

+ Bp{ 1 - p) 

i-p+pexp ycy 


(1 -p) 2 + pexp 

)] 

1 -p + rpexp y^y 



( 22 ) 


This expressions allows us to calculate the bulk particle density in the HD phase. One 
can see that increasing the exit rate (3 lowers the bulk density until the phase boundary 
with the MC phase is reached at p — pmc- This can be illustrated by again considering 
the limiting cases. When motor proteins do not interact with each other (E = 0) Eq. (22) 
yields Phd = 1 — /? and the phase boundary between HD and MC phase can be found 
at ft — 0.5. This is fully consistent with known results for TASEP of non-interacting 
particles 18,19 . For strong repulsions we predict that phd = (1 — /3)/(2 — (3), and 
the phase boundary between HD and MC phases is observed at /3 = 2 — \/2 « 0.59. 
Note that these results slightly differ from calculations for TASEP of dimers because 
of the different exiting rules 30 . For strong attractions the flux through the system is 


vanishing and the high-density is always observed. 

The phase boundary between LD and HD phases can be estimated from the 
condition that at this line the particle currents from both phases become equal, 
Jld = Jhd- It can be shown from Eqs. (16) and ([T7|) that 


— Pld) + (1 — Pld){ 1 — 2 p L o + Pld exp ( -^ff )) 


Phd( 1 - Phd + fP hd exp (-^j ) 


- Phd + Phd exp 
L - Pld + Pld exp 

In this expression, densities phd and p^ are obtained by solving Eqs. (21) and (22), 
respectively, for specific values of E and 6. For E = 0 we find that the LD-HD phase 
boundry is given by (3 = a, and the triple point (where LD, HD and MC phases meet 


together) is found at f3 c = a c = 0.5, as expected for the standard TASEP |h8,19 


.( 23 ) 
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(a) e=o 


° Simulation 
— Theory 


(b) 0=0.25 




Figure 4. Maximal particle currents as a function of the interaction energy for 
different energy splittings: (a) 9 = 0, (b) 9 = 0.25, (c) 9 = 0.75 and (d) 6 = 1. 
In simulations a = /3 = 1 was utilized. 


3. Monte Carlo Simulations and Discussions 

Because of the approximate nature of our method, it is important to test these 
theoretical predictions. It was done in this work by running extensive computer Monte 
Carlo simulations. We utilized the Monte Carlo algorithm known as a random sequential 
update. In our simulations we used a lattice of size L = 1000 to minimize any finite size 
and boundary effects. The particle current and density profiles of molecular motors were 
averaged over 10 8 Monte Carlo steps. To ensure that the system is at the stationary- 
state conditions, the first 20% of events were discarded. We have used a precision of 
0.01 when comparing density profiles to construct accurate phase diagrams. The error 
in calculating the phase boundaries by our method was estimated to be less than 1%. 

In Fig. 3 we compare theoretically calculated phase diagrams with results obtained 
in Monte Carlo computer simulations. One can see that for relatively weak interactions 
theory agrees quite well with computer simulations (Fig. 3c), while for stronger 
interactions (attractive or repulsive) the agreement is mostly qualitative (although still 
better for repulsions): see Figs. 3a, 3b and 3d. Comparing phase behavior at different 
interactions, one can notice that the LD phase is dominating at repulsions, but the 
HD phase is more prevalent for attractions. These observations are consistent with 
expectations that repulsions lead to smaller particle clusters and lower density while 
attractions stimulate the formation of large clusters, which corresponds to higher density. 

However, our method works much better for predicting particle fluxes in the 
MC phase, as illustrated in Fig. 4. The theory correctly describes the fluxes for 
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repulsive interactions for all ranges of parameters (with the exception of the special 
case corresponding to 6 = 0). But for attractions the good agreement is found only 
for small 6. For larger values of the energy splittings (6 > 0.25) there is only a 
qualitative agreement on the overall trends: the fluxes decrease to zero with increasing 
the interaction strength (with the exception of the special case corresponding to 6 = 1). 

These observations suggest that correlations are important for understanding 
dynamic properties of interacting molecular motors. To quantify this effect, we 
investigated a correlation function £7; defined as 


C = (TiTi+ 1) - {n)(r i+ 1). i = 1, ...,L - 1 


(24) 


where two-point and one-point density functions are given by 

(TiT i+ 1 ) = EE TiT i+1 P(Ti,T i+ i ) = Pn, (25) 

Ti T i+1 

(p) = = p- ( 26 ) 

n 


The physical meaning of the correlation function Ci is that it gives a measure of how 
the presence of the particle at site i affects the occupation of the neighboring site i + 1. 
Using the definition together with the normalization condition and Eq. (13), we obtain 
the following analytical expression for C, 

P 2 (l-P) 


C(E) = 


exp 


( —^ - 1 

\k B T 1 


1 + P 


exp 


E 

k B T 


- 1 


(27) 


Note that C is uniform in the bulk of the system. For the case of zero interactions 
(E = 0) it predicts that £7 = 0. This fully agrees with what we know about TASEP 
for noninteracting particles. Here a simple mean-field theory, that completely neglects 
any correlations, works quantitatively well and it correctly describes the majority of 
dynamic properties of the system 118,[l9]. For strong repulsions (E —> —oo) it gives 
C = —p 2 , while for strong attractions (E —>■ oo) we have C = p(l — p). 

The correlation functions predicted in our method and obtained from computer 
simulations are presented in Fig. 5. It is interesting to analyze these data. The physical 
meaning of the correlation function Ci is that it gives a measure of how the presence of 
the particle at site i affects the occupation of the neighboring site i + 1. When there 
are no correlations we have £7 = 0. Negative correlation functions (£7 < 0) indicate that 
there is a less probability to find the particle next to the already occupied site. This is 
the case for repulsive interactions. In contrary, positive values for £7 suggest that the 
presence of the particle at given site enhances the probability to find the particle at the 
neighboring site. It is clear that this situation can be realized for attractive interactions. 
Comparing theoretical predictions with Monte Carlo simulations (Fig. 5) again indicates 
that our theory works very well for repulsive interactions, while for attractions, although 
the trends are correctly picked up, there are deviations. 

The analysis of results presented in Figs. 3, 4 and 5 strongly indicates that 
correlations are important for understanding the mechanisms of interacting molecular 
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(a) 6=0.25 (b) 6=0.5 



Figure 5. Correlations as a function of the interaction energy for: (a) 9 = 0.25, (b) 
9 = 0.5, (c) 9 = 0.75. In simulations a = (3 = 1 was utilized. 


motors. However, it also raises a question of why our theoretical approach, that 
explicitly takes into account some correlations, is able to correctly describe the stationary 
properties only for repulsive interactions and for weak attractions. To answer this 
question we note that the dynamic behavior strongly depends on the sign of the 
interactions. For E < o, the presence of the particle at the site i leads to a lower 
probability of finding another particle at the site i + 1. Then if there is nothing at the 
site i + 1 the occupancy of the site i + 2 will be independent of the fact that there is the 
particle at the site i . These arguments suggest that correlations for repulsive interactions 
are short-range and relatively weak. For E > 0 the situation is very different. Here the 
presence of the particle at the site i stimulates the occupancy of the site i + 1 , and 
consequently the occupation state of the site i + 2 depends on the state of the site i. 
This is consistent with long-range and strong correlations. By construction (see Sec. 
2 . 2 ), our theory accounts only for short-range correlations because the evolution of two- 
site clusters is monitored. This is the main reason why our approach is so successful 
for repulsive interactions, while providing mostly a qualitative description for attractive 
interactions. 

One of the main advantages of our theoretical method is the fact that it can be 
easily extended to account for more realistic features of the motor proteins transport. 
To illustrate this we analyze the effect of varying the splitting coefficient 9 on multi¬ 
particle dynamics of interacting molecular motors. The results are presented in Figs. 4 
and 6 . One can see that dynamics is different for small and large values of the interaction 
splittings. It can be concluded from Eq. (2) that small 9 describe the situation when 
the formation of the particle clusters is weakly affected by the interactions. At the same 
time, the breaking of the particle bonds is strongly influenced by interactions. For 9 « 1 
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the trend is opposite: the particle cluster formation depends strongly on interactions, 
while the bond breaking is not. 

Fig. 4 shows that the particle current (in the MC phase) is generally a non¬ 
monotonic function of interactions. At large repulsions the current saturates, while for 
large attractions the fluxes are going to zero. The maximal particle current is achieved 
for relatively weak repulsive interactions (E ~ —(1 — 2) k#T). The monotonic decrease 
in the particle current is only observed for the special case of 6 = 0. Similar dynamics 
is observed for large 6 > 0.9 (see Fig. 6), but here the most optimal conditions are 
reached now for positive interactions. The case of 6 = 1 is again a special one, and the 
monotonic increase in the current is observed for all range of interactions. 



Figure 6. Maximal particle currents as a function of the interaction energy for 
large energy splittings. Lines are theoretical predictions, symbols are from computer 
simulations. 


In light of these findings, it is important to discuss the effect of intermolecular 
interactions for real motor proteins. These interactions have been measured 
experimentally for kinesins, indicating weak attractions of order of E = (1.6 ± 0.5) 
k#T 10 . Previous theoretical studies suggested that kinesins function at the conditions 


that do support the maximal current, but the analysis was based on the symmetric 
splitting of interactions for transitions rates (6 = 0.5) [l7|. Our new results presented 
in Figs. 4 and 6 indicate that this is probably a reasonable description of multi-particle 
dynamics of kinesins for for most interaction splittings (6 < 0.9). In this case the 
kinesin might operate at the conditions where small changes in interactions lead to 
large modification in the particle dynamics. It has been argued that this might be 


important for maintaining robust cellular transport (5],[8 IT . However, our results (Fig. 
6 ) also suggest another intriguing possibility that the kinesin fluxes might be optimized 
if the splitting affects more the formation of particle clusters (6 > 0.9). The parameter 
6 is a microscopic property that cannot be obtained from our mesoscopic theoretical 
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method. To test this idea it will be important to measure and calculate this quantity 
in more advanced experimental and theoretical investigations. 

4. Summary and Conclusions 

We developed a new theoretical approach to analyze the role of intermolecular 
interactions in the collective dynamics of molecular motors that move along linear 
filaments. Our method is based on utilizing totally asymmetric exclusion processes, 
which have been successfully applied for studying various processes in Chemistry, 
Physics and Biology. It modifies the transition rates by interactions via fundamental 
thermodynamic arguments. A simple theoretical framework, that we call the modified 
cluster mean-held and that takes into account some correlations, is presented and fully 
discussed. It allows us to calculate analytically or numerically exactly all dynamic 
properties of interacting molecular motors. We find that interactions induce correlations 
in the system of collectively moving motor proteins, and the strength of correlations 
depends on the sign of the interactions. It was argued that for repulsions the correlations 
are short-range and relatively weak, while for attractions the range and amplitude of 
correlations are larger. This also leads to different dynamic behavior of interacting 
molecular motors. For repulsions the dynamics is weakly affected by the strength of 
interactions, however for attractions the dynamics is modified much stronger. We also 
investigated the effect of the symmetry of interaction by analyzing splittings between 
different transitions. It was found that when the formation of particle clusters is 
weakly affected by interactions the most optimal fluxes can be realized for weakly 
repulsive interaction. But when breaking of bonds between neighboring particles is 
strongly influenced, the maximal current can be achieved for attractive interactions. 
Furthermore, the importance of these results for kinesin motor proteins has been 
discussed. 

The most important advantage of our method is that it can be easily extended 
to investigate additional more realistic features of molecular motors transport such as 
backward steps, binding and unbinding of motor proteins at all sites, multiple parallel 
pathways and limited resources of motor proteins in the surrounding solution. It will 
be very interesting to generalize our approach to study these phenomena if we want 
to understand better mechanisms of cellular transport processes. However, despite its 
simplicity and the successful application for the interacting molecular motors, it should 
be noted that our method is still approximate and many important features are not 
well described. For example, for a large range of parameters the effect of attractive 
interactions is given only qualitatively. Thus it will be important to test our theory in 
experiments and in more advanced theoretical models. 
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